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Abstract 

A computer code is presented for solving the equations of Hartree-Fock-Bogoliubov (HFB) 
theory by the gradient method, motivated by the need for efficient and robust codes to calculate 
the configurations required by extensions of HFB such as the generator coordinate method. The 
code is organized with a separation between the parts that are specific to the details of the 
Hamiltonian and the parts that are generic to the gradient method. This permits total flexibility 
in choosing the symmetries to be imposed on the HFB solutions. The code solves for both 
even and odd particle number ground states, the choice determined by the input data stream. 
Application is made to the nuclei in the sd-shell using the USDB shell-model Hamiltonian. 
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I. INTRODUCTION 



An important goal of nuclear structure theory is to develop the computational tools 
for a systematic description of nuclei across the chart of the nuclides. There is hardly any 
alternative to self-consistent mean-field (SCMF) for the starting point of a global theory, but 
the SCMF has to be extended by the generator coordinate method (GCM) or other means 
to calculate spectroscopic observables. There is a need for computational tools to carry out 
the SCMF efficiently in the presence of the multiple constraints to be used for the GCM. 
Besides particle number, quantities that may be constrained include moments of the density, 
angular momentum, and in the Hartree-Fock-Bogoliubov (HFB) theory, characteristics of 
the anomalous densities. 

n 

The gradient method described by Ring and Schuck ([1], Section 7.3.3) is very suitable 
for this purpose: it is robust and easily deals with multiple constraints. However, the actual 
computational aspects of the method as applied to HFB have not been well documented 



in the literature. This is in contrast to methods based on diagonalizing the E 



eigenvalue equation. Here there are several codes available in the literature, eg. |^-|8 



FB matrix 



Other, 



less used, methods to solve the HFB equation with multiple constraints can be found in the 
literature; for example the method described in Ref. {q] is close in spirit to the one presented 
here. We note also that the computational issues for using the gradient method in nuclear 
Hartree-Fock theory have been discussed in detail in Ref. That paper also contains 
references to related techniques such as the imaginary time step method. 

Here we will describe an implementation of the gradient algorithm for HFB following 
the iterative method used by Robledo and collaborators [l^. The code presented here, 
hfb_shell, is available as supplementary material to this article (see Appendix). The code 
has separated out the parts that are basic to the gradient method and the parts that are 
specific to the details of the Hamiltonian. As an example, the code here contains a module 
for application to the srf-shell with a shell-model Hamiltonian containing one-body and two- 
body terms. There is a long-term motivation for this application as well. The sd-shell could 
be a good testing ground for the extensions of SCMF such as the GCM and approximations 
derived from GCM. Since one has a Hamiltonian for the sd-shell that describes the structure 
very well, one could test the approximations to introduce correlations, such as projection, the 
random-phase approximation, etc and compare them with the exact results from the Shell 
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12| . As a first step in this 



Model. Preliminary results along this line are discussed in 111 
program, one needs a robust SCMF code that treats shell-model Hamiltonians. Extensions to 
other shell model configuration spaces are straightforward and only limited by the availability 
of computational resources. 

The code described here is more general than earlier published codes in that it can treat 
even or odd systems equally well. The formalism for the extension to odd systems and 



to a statistical density matrix will be presented elsewhere 



13| . We also mention that the 



present code (with a different Hamiltonian module) has already been applied to investigate 



neutron-proton pairing in heavy nuclei 



1^. 



II. SUMMARY OF THE GRADIENT METHOD 

The fundamental numerical problem to be addressed is the minimization of a one- plus 
two-body Hamiltonian under the set of Bogoliubov transformations in a finite-dimensional 
Fock space. We remind the reader of the most essential equations, using the notation of 
Ring and Schuck []]]. The basic variables are the U and V matrices defining the Bogoliubov 
transformation. The main physical variables are the one-body matrices for the density p 
and the anomalous density k, given by 

p = V*V'; K = V*U\ (1) 

The Hamiltonian may be defined in the Fock-space representation as 

H = ei2c\c2 + T H ^1234Ci4c4C3. (2) 
12 ^ 1234 

The expectation value of the Hamiltonian under a Bogoliubov transformation of the vacuum 
is given by 

ifoo = {H) = Tr(£p + \Vp - I A/?*). (3) 

in terms of the fields for the ordinary potential F and the pairing potential A. These are 
defined as 

= X!^1423P34; Ai2 = | X! ^1234't34- (4) 
34 34 

The gradient method makes extensive use of the quasiparticle representation for op- 
erators related to the ordinary and anomalous densities. For a single-particle operator 
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F = ^ij'^l'^j we write 

Y: F,,4c, ^ c^Fc = + /3tFii/3t + 1 [(3F''(3 + ^^F''(3^) . (5) 

where (3, are quasiparticle annihilation and creation operators. The gradients will be 
constructed from the skew-symmetric matrix which for a normal one-body operator is 
given by 

^20 ^ i/ipY* _ v^F^U*. (6) 
The corresponding representation for an operator G of the anomalous density is 

Kc^Gct - cG*c) = 6-°° + + \{(3^G^^(3^ + (3G^^(5) (7) 

The skew-symmetric matrix is given by 

6-20 = U^GU* - V^G*V\ (8) 

Two operators that are particularly useful to characterize the HFB states are the axial 
quadrupole operator Qq and the number fluctuation operator AA^^. We define Qq as 

Qq = 2z^ -x^- y^- (9) 

its expectation value distinguishes spherical and deformed minima. The number fluctuation 
is an indicator of the strength of pairing condensates and is zero in the absence of a conden- 
sate. It depends on the two-body operator A^^, but like the Hamiltonian can be expressed 
in terms of one-body densities. We define it as 

AiV2 ^ ^fj2^ _ ^jY^2 _ Irj.^ (n^^N""^) = 2Tr (p(l - p)) = -2Tr {k*k) . (10) 
The full expansion of the Hamiltonian in the quasiparticle basis is given in Eqs. (E.20- 



E.25) of Here we will mainly need H^^^ given by 

^20 _ ^20 ^ ^20 _ ^\i^y* _ f^t^* _ ^*Y* + f/tA[/*. (11) 

where h = e -|- F. Starting from any HFB configuration U, V one can construct a new 
configuration U', V by the generalized Thouless transformation. The transformation is 
defined by a skew-symmetric matrix Z having the same dimensions as U, V. One often 
assumes that the transformation preserves one or more sjnnmetries such as parity or axial 



rotational symmetry. Then the U, V matrices are block diagonal and Z has the same block 
structure. Otherwise the elements of Z are arbitrary and can be real or complex. The 
transformation is given by 

U' = {U + V*Z*){1 - ZZ*)-^/^ = U + V*Z* + O(Z^) (12) 

V' = {V + U*Z*){1 - zz*)-^i^ = V + \J*Z* + 

The last factor, (1 — ZZ*)"^/^, ensures that the transformed set [/', V satisfies the required 
unitarity conditions for the Bogoliubov transformation. We now ask how the expectation 
value of some bilinear operator Q changes when the Thouless transformation is applied. The 
result is very simple, to linear order in Z: 

= - ^(Tr(Q2°Z*) + h.c.) + (13) 

The same formula applies to the Hamiltonian as well, 

= - ^(Tr(i/^°Z*) + h.c.) + 0{Z^). (14) 

From these formulas it is apparent that the derivative of the expectation value with respect 
to the variables z*^ in Z* is^ 

With a formula for the gradient of the quantity to be minimized, we have many numerical 
tools at our disposal to carry out the minimization. 

It is quite straightforward to introduce constraining fields in the minimization process. 
As seen in Eq. f ll3p the transformation Z will not change the expectation value of Q to 
linear order provided Tr(Q^°Z*) +h.c. = 0. Thus, one can change the configuration without 
affecting the constraint (to linear order) by projecting Z to Zc a.s Zc = Z — XQ^^ with 
A = |(Tr(Q20^*) + h.c.)/Tr(Q20Q2o*)^ ^j^j^ multiple constraints, the projection has the 
form 

Z, = Z-^A„Q^°. (16) 

a 

The parameters are determined by solving the system of linear equations, 

E^-/3Aa = ^(Tr(gJ°Z*) +h.c.) (17) 



^ The derivative is taken with respect to the variables in the skew-symmetric Z*, ie. z*^ — —z*^ and z^j, 



z*j are treated as independent variables. 
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where Map = Tt{Q'^QJ'*). Since we want to minimize the energy, an obvious choice for 
the unprojected Z is the gradient of the Hamiltonian H^^. In this case the constraining 
parameters Aq, are identical to the Lagrange multiphers in the usual HFB equations. We 
will use the notation He for the constrained Hamiltonian 

H, = H-J2^aQa. (18) 
a 

A. Numerical aspects of the minimization 

The most obvious way to apply the gradient method is to take the direction for the change 
from Eq. fll6|17p . and take the length of the step as an adjustable numerical parameter. We 
will call this the fixed gradient (FG) method. It is implemented in the program as 

Z, = r^Hl\ (19) 

Typically the starting V configuration will not satisfy the constraints, and the Z trans- 
formations must also bring the expectation values of the operators to their target values qa- 
The error vector 5qa to be reduced to zero is given by 

^Qa = Ql' - qa. (20) 
We apply Eq. f[T3|) to first order to obtain the desired transformation Z^g, 

Zs, = -Y.M-'p5qaQf. (21) 
With these elements in hand, a new configuration is computed using the transformation 

Z = Z, + Zsq. (22) 

This process is continued until some criterion for convergence is achieved. We shall measure 
the convergence by the norm of the gradient | H^^ \ . This is calculated as 

\Hf\ = {Tr[Hl\Hf)^]f\ (23) 

An example using this method as given is shown in Fig. [H The parameter rj is fixed to some 
value and the iterations are carried out until convergence or some upper limit is reached. 
The required number of iterations varies roughly inversely with 77, up to some point where 
the process is unable to find a minimum in a reasonable number of iterations. 
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FIG. 1: Number of iterations required for convergence using Eq. (jl9p and fixed rj. At the point 
r/ = 0.12 MeV^^ and beyond, the iteration process is unstable. The converged solutions and their 
energies are the same for all values of rj shown in the plot. All values producing converged solutions 
The system is ^^Mg with three constraints, A^, Z, and < Qq >= 10 h/mujQ. The convergence 
criterion is \Hl^\ < 1.0 x 10^^ MeV. See Section IVII Bl for further details. 

There are a number of ways to speed up the iteration process. If the constraints are 
satisfied, the parameter rj can be increased considerably. Fig. [2] shows the change in H^^ from 
one iteration cycle as a function of 1] using to update. For small values of rj, the change 
in constrained energy is given by the Taylor expansion Eq. (IT^ . AH^^ ^ —rjTr 
This function is shown as the straight line in the Figure. The actual change is shown by 
the black circles. One sees that rj could be doubled or tripled from the maximum value 
permitted in Fig. [TJ However, the constraints and other aspects of the new U, V become 
degraded so that such steps are not permissible for many iterations [2]. Still, one can take 
advantage of the possible improvement by choosing t] at each iteration taking account of the 
relevant information from the previous iteration. This can be extracted from the ratio 

A/700 

r = (24) 

which is close to one for too-small rj values and close to | at the value corresponding to the 
steepest-descent minimum. We call such methods variable gradient. We note that updates 
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FIG. 2: Single-step energy change as a function of rj in Eq. (jl9p . The configuration that was 
updated is the 10th iteration step of the system in Fig. 1. 

with Zsq alone are relatively quick because there is no need to evaluation matrix elements 
of the Hamiltonian. These considerations are implemented in the code of Ref. 10| by 
interspersing cycles of iteration by Zsg alone among the cycles with updates by Eq. fl22|) . 

Another way to improve the efficiency of the iteration process is to divide the elements 
of iff by preconditioning factors Pij, 



(25) 



The choice of the preconditioner is motivated by Newton's method to find zeros of a 
function (here iff) based on knowledge of its derivative. This could be accessible from 
the second-order term in Eq. f|T^ . but unfortunately it cannot be easily computed as 
it involves the HFB stability matrix. However a reasonable approximation to it can be 
obtained from H^^, the one-quasiparticle Hamiltonian that, when in diagonal form, is the 
dominant component of the diagonal of the stability matrix. One first transforms U,V to a 
basis that diagonalizes H^^. Call the eigenvalues of the matrix Ei and the transformation 
to diagonalize it C. The U, V are transformed to U', V in the diagonal quasiparticle basis 
by 

U' = UC] V = V'C (26) 
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In the new basis the preconditioner is given by 

Pij = max^Ei + Ej, Emin) (27) 

where Emin is a numerical parameter of the order of 1-2 MeV. The main effect of the 
preconditioner is to damp away those components of the gradient with high curvatures (i.e. 
second derivatives) which correspond to two-quasiparticle excitations with large excitation 
energies. This is very important for Hamiltonians that have a large range of single-particle 
energies, such as the ones derived from commonly used nuclear energy density functionals 
such as Skyrme and Gogny. 

In Table I we show the number of iterations required to reach convergence for a case 
calculated in Table II, to be described below. We see that there is a gain of more than 



Method 


^ Vmin Vniax 


Iconv 


fixed gradient 


0.10 MeV"^ 


140 


variable gradient 


0.08 MeV"^ 0.3 MeV"^ 


65 


fixed pr. 


0.7 


72 


variable pr. 


0.7 2.0 


34 



TABLE I: Number of iterations to convergence Iconv with various treatments of the update. Eq. 
([T9|) with fixed and variable gradients is used for the top two lines and the preconditioned gradients 
Eq. (|25p are used for the lower two lines. The system is ^^Ne as calculated in the top first entry 
in Table II. 



a factor of 3 between the naive steepest descent and the preconditioned gradient with a 
variable rj. Similar ideas have been used in a HF context in {2], [isl with similar speedups. 

III. ODD-A NUCLEI 

As discussed by Ring and Schuck[l|, each U^V set can be characterized by its number 
parity, either even or odd. This means that when the wave function is constructed and states 
of definite particle number are projected out, the nonzero components will have either all 
even or all odd particle number. Another important fact is that the generalized Thouless 
transformation does not change the number parity of the Bogoliubov transformation. Thus, 
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if we start from aU,V set of odd number parity, the final converged configuration will only 
have components of odd nucleon number. 

In fact, in the matrix-diagonalization method of solving the HFB equations, the higher 
energy of the odd-A configurations requires some modification to the Hamiltonian or to the 
iteration process. A common solution is to add additional constraining fields so the that 
odd-A system has lower energy JigI. 'v\ . Typically the external field to be added breaks time 
reversal symmetry in some way. But then one can no longer assert that a true minimum has 
been found, because the extra constraints can affect the configuration. The gradient method 
does not have this shortcoming. If the space of odd-number parity Bogoliubov transforma- 
tions is adequately sampled, it will find the global minimum of the odd-A configurations. 
Moreover, with the gradient method one does not need to modify the computer code to treat 
odd-A systems. Only the initial U, V set is different for the two cases. 

We note the H^^ has negative quasiparticle eigenenergies in the odd number-parity space, 
assuming that the true minimum of the HFB functional is an even number-parity configu- 
ration. 



IV. OTHER SPECIAL CASES 



The variational minimum might not be directly reachable by the generalized Thouless 
transformation, but it always is a limit of a succession of transformations. This is the case 
if the condensate vanishes at the minimum while the starting configuration has a finite 
condensate. This does not cause any practical difficulties except for reducing the rate of 
convergence. Still, in such cases it is more direct to start with a. U,V configuration of the 
pure Hartree-Fock form. It is not possible to use the gradient method in the other direction, 
to go to a minimum having a finite condensate from a starting U, V of Hartree-Fock form, 
as explained below. 



V. IMPOSED SYMMETRIES 



The U, V matrices have a dimension of the size of the Fock space of nucleon orbitals and 
in principle can be dense matrices. However, one often imposes symmetries on the wave 
function by assuming that the U, V have a block structure with all elements zero outside 
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the blocks. For example, most codes assume separate blocks for neutrons and protons. This 
is well-justified when there is a significant difference in neutron and proton numbers but 
in general it is better to allow them to mix. Other quantum numbers that are commonly 
imposed on the orbital wave functions are parity and axial symmetry. There are only a 
few exceptional nuclei that have HFB ground states breaking these symmetries. For the 
parity, there are the Ra nuclei and Th nuclei. Concerning axial symmetry, a global study 



found only three cases of nonaxial HFB 



of even-even nuclei with the Gogny functional 
minima among 1712 nuclei. 

The number of orthogonal minima that can be easily calculated in the gradient method 
depends on the assumed block structure. In the even number-parity space there is just one 
global minimum. But in the odd number-parity space the number parity of each block is 
conserved in the iteration process, so there will be one state for each block. For example, 
states of different i^-quantum number may be calculated by imposing a block structure 
that imposes axial symmetry. Thus for odd-A nuclei, the quasiparticle can be in any of the 
i^'-blocks, giving a spectrum of states with K specified by the block. 

A more subtle form of possible imposed symmetries is those contained in the starting 
U, V configuration. The energy H^^ is essentially a quadratic function of symmetry-breaking 
densities because the products of densities in the functional must respect the symmetries 
of the Hamiltonian. If these components are zero in the initial configuration, the energy 
is stationary at that point and there is no gradient to generate nonzero field values. The 
typical cases are quadrupole deformation in the ordinary density and any form of anomalous 
densities. Fortunately, it is very easy to avoid unwanted symmetries in the starting U, V as 
discussed below. 



VI. THE CODE HFB_SHELL 

The code hfb_shell presented in this paper is described in more detail in the Appendix. 
The main point we want emphasize about the code is that it is organized in modules that 
separate out the functions that are independent of the Hamiltonian from those that are 
specific to it. Also, the block structure is specified only by the code input, and can easily be 
changed. The examples we show are for the sci-shell using the USDB Hamiltonian jl^. Since 
that Hamiltonian is specified by the fitted numerical values of the 3 single-particle energies 
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and the 63 JT-coupled two-particle interaction energies, it does not have any symmetries 
beyond those demanded by the physics. In particular, the HFB fields obtained with it 
should provide a reahstic description of aspects such as the time-odd fields, that are difficult 
to assess with the commonly used energy functionals such as those in the Skyrme family. 

A. Application to the sd-shell 

The sd shell- model space has a dimension of 24 and the principal matrices U,V,Z,... have 
the same dimension. In the application presented here, we assume axial symmetry which 
splits the matrices in blocks of dimension 12, 8 and 4 for m-quantum numbers ±|, ±|, 
and ±1 respectively. Neutron and proton orbitals are in the same blocks, so the basis is 
sufficiently general to exhibit neutron-proton pairing, if that is energetically favorable. We 
also assume that the matrices are real. 

We often start with a. U,V configuration of canonical form, namely U diagonal, Uij — 
u5ij. The nonzero entries of the V are all equal to ±v — ±-\/l — u'^, and are in positions 
corresponding to pairing in the neutron-neutron channel and the proton-proton channel. 
We arbitrarily take u = 0.8 and v = 0.6 for the starting configuration Uq^Vq. This may be 
modified in a number of ways before it is used as a starting configuration in the gradient 
minimization. When calculating a nucleus for which or Z is zero or 12, it is more efficient 
to use U, V matrices that have those orbitals empty or completed filled in the starting 
configuration. This is carried out by changing u, v to zero or one for the appropriate orbitals. 
The particle number of that species is then fixed and is not constrained in the gradient search. 

For odd-number parity configurations, the U, V is changed in the usual way by inter- 
changing a column in the U matrix with the corresponding column in V. The space that 
will be searched in the gradient method then depends on the block where the interchange 
was made. In principle it does not depend on which column of the block was changed. 
However, there is some subtlety is making use of this independence which will be discussed 
below. 

We may also apply a random Z transformation to the starting configurations. Since all 

the entries in the upper triangle of the Z matrix are independent, we can populate them 
with random numbers. This seems to be a good way to break unwanted symmetries in the 
starting configuration that would be preserved by the gradient update. We denote by Ur,Vr 
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the configuration generated from f/o, Vq by a randomly generated Z. 

In principle one could also start from the [/, V configuration of the vacuum: f/ = 1, K = 0. 
We have tried this and found, as might be expected, that the proportion of false minima is 
larger than is obtained with f/g, Vq. 



VII. THREE EXAMPLES 

In this section we will describe the HFB calculations for three nuclei, ^^Mg, ^^Mg, and 
^^Ne. The first one is typical of a spherical nucleus that exhibits identical-particle pairing. 
The second is a well-deformed nucleus. The third illustrates the method for an odd-A 
system. 

For calculating matrix elements of the quadrupole operator Qq, we will treat the single- 
particle wave functions as harmonic oscillator functions of frequency lo^^ and report the 
quadrupole moments in units of Ti/muQ. 
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Mg 

The nucleus ^^Mg {{N,Z) = (12,4) in the S(i-shell) behaves as expected of a semimagic 
nucleus in HFB. Please note that we do not include in our configuration space the /7/2 



2Q|, 



2l|. We 



intruder shell required to explain the deformation properties of this nucleus 
calculate the HFB ground state in two ways, illustrating the role of the starting configuration. 
The first is to use a randomized Ur, K configuration and constraining the particle numbers 
to the above values. Another way is to start with a prolate configuration similar to Uo,Vq 
for the protons and with all the neutron orbitals filled. In that case, only the proton number 
is constrained. Both iteration sets converge to the same minimum, a spherical configuration 
having a strong proton pairing condensate. The output characteristics are Ehfb = —135.641 
MeV, Qq = 0.00 and AZ^ = 2.93. The zero value for shows that the configuration is 
spherical, and the nonzero value for AZ'^ shows that protons are in a condensate. Next we 
calculate the condensation energy, defined as the difference between Ehfb and the Hartree- 
Fock minimum Ehf- The easiest way to find the HF minimum is to repeat the calculation 
with an additional constraint that forces the condensate to zero. This is done by adding a 
G-type operator that is sensitive to the presence of a condensate. Carrying this out, we find 
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a minimum at Ehf = —134.460 MeV and = 5.08. The extracted correlation energy is 
Ehf — Ehfb = 1-18 MeV, which is much smaller than what one would obtain with schematic 
Hamiltonians fitted to pairing gap. It is also interesting to extract the quasiparticle energies, 
since they provide the BCS measure of the odd-even mass differences. These are obtained 
by diagonalizing H^^. The results for the HFB ground state range from 1.5 to 9 MeV, with 
the lowest giving the BCS estimate of the pairing gap. 

B. 24]yig 

The next nucleus we consider, ^"^Mg with = 4 and Z = 4, is strongly deformed in 
the HFB ground state. We find that the converged minimum has a quadrupole moment 
(Qq) = 12.8, close to the maximum allowed in the space. More surprisingly, the pairing 
condensate vanishes at the HFB convergence. We now make a set of constrained calculations 
to display the energy as a function of quadrupole moment. The starting configuration is 
generated by applying a random transformation to Uq, Vq. The gradient code carries out the 
iterations with the constraints N = A, Z = 4, and the chosen value of Q. The convergence 
of the constraints to their target values is very rapid, using the update in Eq. fl2T]) . This is 
illustrated in Fig. |3l showing the deviation from the target values as a function of iteration 
number in one of the cases {Q = 10). On the other hand, the convergence to the minimum 
of the HFB energy can be slow, using a fixed- update with Eq. (IT^ . The calculations were 
carried out setting the convergence criterion \H^^\ < 0.01 MeV. Fig. |4]shows the number of 
iterations required to reach convergence for the various deformations. They range from 40 
to 250. In a number of cases, the iterations seem to be approaching convergence, but the 
system is actually in a long valley, and eventually a lower minimum is found. It may also 
happen that the gradient method finds a local minimum that is not the global one. Perhaps 
10% of the runs end at a false minimum. This can often be recognized when carrying 
constrained calculations for a range of constraint values, as it gives rise to discontinuities in 
the energy curves. The only systematic way we have to deal with the false minima is to run 
the searches with different randomly generated starting configurations, and select the case 
that gives the lowest energy. The resulting deformation plot combining two runs is shown in 
Fig. O The global minimum is at a large prolate deformation as mentioned earlier. There is 
also a secondary minimum at a large oblate deformation. For all deformations, the ordinary 



14 



10 



< 



0.1 



0.01 




0.001 ^ ' ' ' ^ 

12 3 4 

Iteration number 

FIG. 3: Error in constrained quantities as a function of iteration number for the ?7 = 0.1 run of 
the ^^Mg iterations in Fig. 1. Quantities constrained are: N, open circles; Z, filled squares; and 
Qq, filled circles. 
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FIG. 4: Number of iterations required to convergence for the calculated configurations on the 
deformation energy curve Fig. [5l 
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FIG. 5: HFB energies as a function of deformation, using the Qq quadrupole constraint. The 
nucleus is ^'^Mg, = Z = 4 in the sd-shell. 



neutron-neutron and proton-proton pairing condensates are small or vanish. 



C. 2iNe 

The next nucleus we discuss, ^""^Ne with {N,Z)sd = (3,2), illustrates how the gradient 
method makes use of the conserved number parity to find the minimum of odd-A systems. 
We start with the Uq, Vq configuration, and convert it to an odd-number parity configuration 
by exchanging two columns in the m = ±| block. There are 6 possible columns with m = +| 
that can be exchanged. The results for the converged energies are shown in the top row of 
Table Ull All of the neutron exchanges give the same final energy, —40.837 MeV. However, 
the energy is different for proton exchanges. The reason is that the starting configurations do 
not mix neutrons and protons, and for reasons discussed earlier the corresponding gradients 
are zero. This unwanted symmetry can be broken by making a random transformation of 
the initial configuration. The results are shown in the second row. Now all the energies are 
equal, showing that the minimum can be accessed from any column exchange. Interestingly, 
the energy is lower than in the previous set of minimizations. This shows that there is a 
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significant neutron-proton mixing in the condensate for ^^Ne. 



u,v 


"5/2,1/2 "3/2,1/2 


„n JP 
*l/2,l/2 "5/2,1/2 


d^ 

"3/2,1/2 '^1/2, 1/2 


Uo,Vo 

Ur,Vr 


-40.837 -40.837 
-41.715 -41.715 


-40.837 -40.215 
-41.715 -41.715 


-40.176 -40.176 
-41.715 -41.715 



TABLE II: HFB energies of Ne, with different starting configurations. For the top row, the 
starting configuration is Uo,Vo with the indicated column in the m = ±^ block interchanged. The 
second row starts from a randomized configuration Ur,Vr as discussed in Sect. IVI A[ 
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Appendix: explanation of the code 

The code hf b_shell that accompanies this article implements the gradient method dis- 
cussed in the text^. The code is written in Python and requires the Python numerical library 
numpy to run (see 



221 ] and accompanying papers for a description of Python in a scientific 
environment). The main program is the file hfb.py. It first carries out the initialization 
using information from the primary input data file that in turn contains links to other 
needed data files. There are three of these, one for the Hamiltonian parameters, one for the 
correspondence between orbitals and rows of the U, V matrices include the assumed block 
structure, and one for the input U, V configuration. The input data format is explained in 
the readme . txt of the code distribution. 

Following initialization, program enters the iteration loop, calling the various functions 
used to carry out the iteration. The loop terminates when either a maximum number 
of iterations itmax is reached or the convergence parameter \H^^\ go below a set value 
converge. 



^ The code may be downloaded from http : //www. phys .Washington, edu/users/bertsch/hfb-shell . 21 .tar 

until it has been published in a journal repository. 
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The function calls that are specific to the sd-shell application are collected in the module 
sd_specif ic .py. The tasks carried out by these functions include: 

• initialization of matrix sizes and block structures 

• setting up the matrices representing single-particle operators in the shell-model basis. 

• calculation of the fields F, A from the densities p, k. This function makes use of a 
table of interaction matrix elements Vijki that are read in from a file. The present 
distribution of the code only provides the Hamiltonian data for the USDB interaction 

m- 

The functions that are generic to the gradient method are collected in the module 
hf b_utilities .py. Many of these functions are defined by equations in the text; the 
correspondence is given in Table III. 



Function call 


Equation in text 


rho_kappa 


dlD 


F20 


m 


G20 


m 


H20 


m 


HOO 


m 


Ztransf orm 


m 



TABLE III: Python functions in hfb_utilities .py corresponding to equations in the text. 



The output of hf b . py reports the expectation values of the Hamiltonian and the single- 
particle operators A^, Z and Qq at each iteration step, together with the convergence pa- 
rameter \H^^\. After the final iteration, the values are reported for the expectation values 
of constraining parameters and the number fiuctuations AN"^, AZ"^. The final U, V con- 
figuration is written to the file uv.out. Thus additional iterations can be performed simply 
by specifying uv.out as the new input file. 

In addition, there is a set of functions collected in the module hf b_tools .py. These are 
useful for making input U, V configurations and for analyzing the output U, V configuration, 
but are not needed to run hf b . py. For example, a randomizing transformation can be applied 
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to a, U,V configuration by the function randomize. Another useful function is canonical, 
used to extract the eigenvalues of the p operator needed for the canonical representation. 
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